%% master_counterfactuals_CES_approx_countryloop.m


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Author: Produced for JdG,AL,IM by Christopher Evans at UPF 
%
% Date: February 2019
%
% Purpose: This file runs the counterfactual exercises after the
% main_baseline.m has been run. Change the shocks to experiment with the
% counterfactuals.
% 
% Based on one file: Step_4_Algorithm_counterfactual.m, which uses new
% baseline input data, MAT/Step_4_setup_`alphaT''rhoT'_`shockT'.mat
% as initial level variables, and then shocks the system and computes hat
% variables for wages, prices, exports, etc.
%
% NEW: for CES setup. Have streamline file that has to be called up and
%                     various variables that have to be used.
% Output: MAT/Step_4_results_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat,
%         MAT/Step_5_results_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat
%         MAT/Elasticity_output_`alphaT''rhoT'`lambdaT'`etaT'_`shockT'.mat
%
% Last Updated: 07/09/2019 by JDG/IM
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

global alphaT rhoT lambdaT etaT shockT rho sigma sigmaT lambda eta ...
    tol maxit Pi_hat_n D_hat_n Xi_hat_mnjf Xi_hat_mnj a_hat_f a_hat ...
    varphi sL_n sPi_n sD_n N J FI_J france oneminus_alpha_vector ...
    labour_share_PWT countrysecD firmsecD_sorted start start_sorted ...
    sum_by_country_dummy ISO alpha_WIOD_j alpha_WIOD_francej sector_algorithm 

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% alphaT: 
%
% 1 = PWT labor share for all countries and Labor share for French firms 
%     are at the COUNTRY level.
%
% 2 = PWT labor share for all countries and WIOD scaled labor share for 
%     French firms. Labor share for French firms are at the FIRM level.
%
% 3 = WIOD labor share for all SECTORS with the labor share at the SECTOR 
%     level for French firms.
%
% 4 = WIOD labor share for all SECTORS with the labor share at the FIRM 
%     level for French firms 
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

alphaT = 4;

%% Loading in the previous data to get some useful variables

fname = strcat('MAT/Step_2_firm_data_alpha',num2str(alphaT),'.mat');
load(fname,'FI_J', 'firmsecD_sorted', 'start_sorted', 'alphaT');
load('MAT/WIOD.mat','countrysecD','france','J','N','start','sum_by_country_dummy');

%% SET PARAMETERS VALUES

% Frisch elasticity for GHH preferences

varphi=3;

% share of labour in total final consumption (need to make assumption or collect data)

sL_n = ones(1,N)*0.8; 

% share of profits in total final consumption (should have this data)

sPi_n= ones(1,N)*0.15; 

% share of trade deficit in total final consumption (should have this data)

sD_n= ones(1,N)*0.05; 



%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% rhoT:
%
% Value for rho. Will apply to J number of sectors to create a vector of
% homogeneous rhos across sectors.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

rhoT = '3';

rho = 3*ones(J,1);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Value for sigma
%
% This is the elasticity of substitution across sector consumption. Given
% we do not adjust this anymore after running C-D production we set to the
% value we think is most reasonable: 1.5. There is therefore no longer a
% "type" variable.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

sigma = 1.5*ones(J,1);
sigmaT = '1.5';

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% lambdaT:
%
% Value for lambda. We might experiment on this so have a type.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

lambdaT = '1.000001';

lambda = 1.000001;

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% etaT:
%
% Value for eta. We might experiment on this so have a type.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

etaT = '1.000001';

eta = 1.000001;

%% Baseline: no shocks

% Pi_hat: Change in profits per country
    
Pi_hat_n=ones(1,N);
    
% D_hat: Change in trade deficit
    
D_hat_n = ones(1,N);
    
% Taste shock - for France only
    
Xi_hat_mnjf=ones(FI_J,N);
    
%Taste shock - for sectors only (countries except France)
    
Xi_hat_mnj= ones(N*J,N);
    
% Productivity shock, a is firm specific - France only
    
a_hat_f= ones(FI_J,1);
    
% Productivity shock, a is for sectors - countries except France
    
a_hat = ones(J*N,1);

ISO = {'AUS', 'AUT', 'BEL', 'BGR', 'BRA', 'CAN', 'CHN' ,'CYP', ...
           'CZE', 'DEU', 'DNK', 'ESP' , 'EST', 'FIN', 'FRA', 'GBR', ...
           'GRC', 'HUN', 'IDN', 'IND', 'IRL', 'ITA', 'JPN', 'KOR', ...
           'LTU' ,'LVA', 'MEX','MLT', 'NLD', 'POL', 'PRT', 'ROU','RUS', ...
           'SVK', 'SVN', 'SWE', 'TUR', 'TWN', 'USA','ROW'};
    

for qq = 1:1
    %1:4 'BEL', 'DEU',
    %'AUS' 
    for cty = {'AUT' 'BEL' 'BGR' 'BRA' 'CAN' 'CHN' 'CYP' ...
           'CZE' 'DEU' 'DNK' 'ESP'  'EST' 'FIN' 'GBR' ...
           'GRC' 'HUN' 'IDN' 'IND' 'IRL' 'ITA' 'JPN' 'KOR' ...
           'LTU' 'LVA' 'MEX' 'MLT' 'NLD' 'POL' 'PRT' 'ROU' 'RUS' ...
           'SVK' 'SVN' 'SWE' 'TUR' 'TWN'  'USA' 'ROW'};
    idx = find(ismember(ISO, cty));
%% Feeding in SHOCKS
    
    if qq == 1
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Shock 1:
    %
    % A +10% productivity shock to all firms for cty. 
    % UNIQUE IDENTIFIER (shockT) keeps track of counterfactual exercise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    shockT = strcat(cty,'prod_p10');
    shockT=shockT{1};
    a_hat = ones(J*N,1);
    a_hat(start(idx):start(idx+1)-1)=1/1.1;
    Xi_hat_mnjf=ones(FI_J,N);
    
    elseif qq == 2 
      %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    % Shock 2:
    %
    % A +10% preference demand to all firms from cty. 
    % UNIQUE IDENTIFIER (shockT) keeps track of counterfactual exercise
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    shockT = strcat(cty,'pref_p10');
    shockT=shockT{1};
    Xi_hat_mnjf=ones(FI_J,N);
    Xi_hat_mnjf(:,idx) = 1.1^(str2num(rhoT)-1);
    a_hat = ones(J*N,1);
    
    end
  
    %% Running the Algorithm
    
    disp('BEGINNING COUNTERFACTUAL EXERCISE')
    
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%   
    % Step_4_Algorithm: same algorithm as Step_4_Algorithm_baseline.m, but now
    % we are able to feed in the model the calculated t+1 varialbes as t,
    % starting at the new equilibrium with no wedges
      eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_setup_alpha',num2str(alphaT),...
      '_rho',rhoT,'_approx.mat')]); 
       Step_4_Algorithm_counterfactual_CES_approx
% %     
%     % Save the counterfactual output so it can be used for post analysis
% 
    outputtable = [' X_mnj X_hat_mnj pi_mnjf1 pi_mnj1 ' ...
           'pi_c_mnj1 pi_c_nj1 pi_M1 pi_l1 pi_M_f1 pi_l_f1 '...
           'PC_hat_n1 PC_n P_hat_n P_hat_nj P_hat_mnj w_hat ' ...
           'a_hat a_hat_f Xi_hat_mnj Xi_hat_mnjf '...
           'final_goods1 intermediates1'] ;
% % %   
    eval([strcat('save MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_results_alpha',num2str(alphaT),...
       '_rho',rhoT,'_',shockT,'_approx.mat',outputtable,' -v7.3;')]); 
% 
% %     
    display('STEP 4 COUNTERFACTUAL COMPLETED SUCCESFULLY')
%     
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    
    % Step_5_counterfactual_output: Allows the user to easily analyse the
    %      counterfactual exercises produced before. Runs the CounterOutput.m
    %      function after specifying alpha type (labor share), gamma
    %      calibration and the shock scheme to be analysed.               %% Load data

            
            % Time 0 data

             eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,...
     '_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_setup_alpha',num2str(alphaT),...
      '_rho',rhoT,'_approx.mat')]);
% 
%             % Time 1 data
% 
             eval([strcat('load MAT/rho',rhoT,'eta',etaT,'lambda',lambdaT,...
      '_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_4_results_alpha',num2str(alphaT),...
      '_rho',rhoT,'_',shockT,'_approx.mat')]);
%  
%     
      Step_5_counterfactual_output_CES_approx
%     
%     % Save the counterfactual output so it can be used for post analysis
% 
 outputtable = [' VA_hat_n VA_hat_mj VA_hat_mnj VA_hat_fj VA_hat_fnj '...
           'RGDP_def_hat_n RGDP_def_hat_fj RGDP_cpi_hat_n RGDP_cpi_hat_fj '...
           'RGDP_doubledef_hat_n RGDP_doubledef_hat_fj '...
           'real_wage_cpi_hat real_wage_def_hat real_wage_doubledef_hat '...
           'real_wage_income_cpi_hat_n real_wage_income_def_hat_n real_wage_income_doubledef_hat_n '...
           'TFP_hat_def_n TFP_hat_cpi_n TFP_hat_doubledef_n TFP_hat_cpi_france_j '...
           'imports_over_GDP_all imports_over_GDP_france_sector '...
           'exports_over_GDP_all exports_over_GDP_france_sector '...
           'GDP_deflator_n DoubleDeflation_deflator_n P_hat_n VA_fj_0 VA_mj_0'] ;
       'TFP_hat_def_france_j TFP_hat_doubledef_france_j '...
% %  
      eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
          '_rho',rhoT,'_',shockT,'_approx.mat',outputtable,' -v7.3;')]); 
%   
% %     
      display('STEP 5 COMPLETED SUCCESFULLY')
%     
    %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
    
    % Step_6_Elastisticity_output: perform covariance decomposition at firm and
    %        sector level

    % Load files: only needed if run separately from Step 5

    eval([strcat('load MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_5_output_alpha',num2str(alphaT),...
         '_rho',rhoT,'_',shockT,'_approx.mat;')]); 

    % Shock size
    
   delta = 0.1;
    
   Step_6_Elasticity_output_CES

             outputtable = [' e_Y_FRf_cpi m_e_f_cpi cov_e_f_cpi check_f_cpi '...
           'e_Y_FRs_cpi m_e_mjFRs_cpi cov_e_mjFRs_cpi check_mjFRs_cpi '...
           'e_Y_s_cpi m_e_mj_cpi cov_e_mj_cpi check_mj_cpi '...
           'e_Y_FRf_def m_e_f_def cov_e_f_def check_f_def '...
           'e_Y_FRs_def m_e_mjFRs_def cov_e_mjFRs_def check_mjFRs_def '...
           'e_Y_s_def m_e_mj_def cov_e_mj_def check_mj_def '...
           'e_Y_FRf_doubledef m_e_f_doubledef cov_e_f_doubledef check_f_doubledef '...
           'e_Y_FRs_doubledef m_e_mjFRs_doubledef cov_e_mjFRs_doubledef check_mjFRs_doubledef '...
           'e_Y_s_doubledef m_e_mj_doubledef cov_e_mj_doubledef check_mj_doubledef'] ;
%  

     eval([strcat('save MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT),'/Step_6_output_alpha',num2str(alphaT),...
         '_rho',rhoT,'_',shockT,'_approx.mat',outputtable,' -v7.3;')]); 

     display('STEP 6 COMPLETED SUCCESFULLY')
    
    disp('COUNTERFACTUAL EXERCISE COMPLETED SUCCESFULLY')
    end
end
 
% Compare multi-country and single-country elasticities
%fname = strcat('MAT/rho',num2str(rhoT),'eta',num2str(etaT),'lambda',num2str(lambdaT),'_psi',num2str(varphi),'_sigma',num2str(sigmaT)/Step_5_output_alpha4_rho3_WORLDprod_p10_approx.mat');   
%load(fname,'VA_fj_0','VA_hat_fj');

World_VA_hat_fj=VA_hat_fj-1;

CTY = {'AUS', 'AUT', 'BEL', 'BGR', 'BRA', 'CAN', 'CHN' ,'CYP', ...
'CZE', 'DEU', 'DNK', 'ESP' , 'EST', 'FIN', 'FRA', 'GBR', ...
'GRC', 'HUN', 'IDN', 'IND', 'IRL', 'ITA', 'JPN', 'KOR', ...
'LTU' ,'LVA', 'MEX','MLT', 'NLD', 'POL', 'PRT', 'ROU','RUS', ...
'SVK', 'SVN', 'SWE', 'TUR', 'TWN', 'USA','ROW'};

ByCountry_VA_hat_fj=zeros(FI_J,1);

for i = 1:size(CTY,2)
   cty = string(CTY(i));
   if cty~='FRA' & cty~='ROW'
   % Get required counterfactual data
   
   fname = strcat('MAT/rho3eta1.000001lambda1.000001_psi3_sigma1.5/Step_5_output_alpha4_rho3_',cty,'prod_p10_approx.mat');   
   load(fname,'VA_hat_fj');   
   ByCountry_VA_hat_fj=ByCountry_VA_hat_fj+(VA_hat_fj-1);
   end
end

scatter(World_VA_hat_fj,ByCountry_VA_hat_fj)
hold on
plot(World_VA_hat_fj,World_VA_hat_fj)
ylabel({'Country-specific shocks'});
xlabel({'Multi-country Shock'});
hold off